/*
 * Assembler.cpp
 *
 *  Created on: 7 Sep 2011
 *      Author: Allan
 */

#include "Assembler.h"

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "Geochemistry/GeoDatabase.h"

#include "Core/Phase.h"
#include "Core/Multiphase.h"
#include "Core/Reaction.h"
#include "Core/ReactionSystem.h"

#include "Utils/Algorithms.h"
#include "Utils/SpeciesUtils.h"

const vector<string> AssembleSpecies(const vector<string>& components)
{
	vector<string> primarySpecies = AssemblePrimarySpecies(components);
	
	vector<string> secondarySpecies = AssembleSecondarySpecies(components);
	
	return Union(Union(primarySpecies, secondarySpecies), components);
}

const vector<string> AssembleAqueousSpecies(const vector<string>& components)
{
	return FilterSpecies(AssembleSpecies(components), {"aqueous"});
}

const vector<string> AssembleGaseousSpecies(const vector<string>& components)
{
	return FilterSpecies(AssembleSpecies(components), {"gaseous"});
}

const vector<string> AssembleMineralSpecies(const vector<string>& components)
{
	return FilterSpecies(AssembleSpecies(components), {"mineral"});
}

const vector<string> AssemblePrimarySpecies(const vector<string>& components)
{
	// The set containing only unique basis species
	set<string> primarySpecies;
	
	// Make sure the hydrogen ion H[+] is present in the set of primary species in case H2O is listed in (components)
	if(IndexOf("H2O", components) < components.size())
		primarySpecies.insert("H[+]");
	
	// Loop over all the provided (components)
	BOOST_FOREACH(const string& s, components)
	{
		// Check if the current known species is a basis species in the database
		if(GeoDatabase::GetInstance().GetSpeciesInfo(s).HasTag("basis"))
			primarySpecies.insert(s);
		else
		{
			// Get the info data of the reaction whose reactant is the current known species
			const ReactionInfo reactionInfo = GeoDatabase::GetInstance().GetReactionInfo(s);
		
			// Add all the products in the reaction to the set of basis species being created (all these products are basis species in the database)
			primarySpecies.insert(reactionInfo.products.begin(), reactionInfo.products.end());
		}
	}
	
	// Convert the set of basis species to a vector container and return it
	return vector<string>(primarySpecies.begin(), primarySpecies.end());
}

const vector<string> AssembleSecondarySpecies(const vector<string>& components)
{
	const vector<string> primarySpecies = AssemblePrimarySpecies(components);
	
	// Get a reference to the list of reactions in the database
	const GeoDatabase::ReactionInfoList& reactionInfoList = GeoDatabase::GetInstance().GetReactionInfoList();
	
	// The vector of secondary species that can be built from the set of basis species
	vector<string> secondarySpecies;
	
	// Loop over all the reactions in the database
	BOOST_FOREACH(const ReactionInfo& r, reactionInfoList)
	{
		// An auxiliary boolean variable to determine if all the products of the current reaction belong to the set of basis species
		bool belong = true;
		
		//  Check if the all the products of the current reaction belong to the set of basis species (basisSpecies)
		BOOST_FOREACH(const string& product, r.products)
			if(!IsContained(product, primarySpecies))
				belong = false;
		
		// Insert the reactant (a secondary species) of the current reaction into the vector of secondary species
		if(belong)
			secondarySpecies.push_back(r.reactant);
	}
	
	return secondarySpecies;
}

const Reaction AssembleReaction(const string& reactant, const Multiphase& multiphase)
{
	// Get the information data of the reaction whose reactant is (reactant)
	const ReactionInfo rinfo = GeoDatabase::GetInstance().GetReactionInfo(reactant);
	
	// The species that participate in the reaction
	vector<string> species = rinfo.products; 
	species.push_back(rinfo.reactant);
	
	// The stoichiometries of the species in the reaction
	vector<double> stoichiometries = rinfo.stoichiometries; 
	stoichiometries.push_back(-1.0);
	
	// Create the reaction equation
	ReactionEquation equation(species, stoichiometries);
	
	// Create the equilibrium constant of the reaction
	EquilibriumConstant K(rinfo.logKi, rinfo.Ti);
	
	// Create the reaction
	Reaction reaction(equation, K, multiphase);
	
	return reaction;
}
